Internal dissipation of a polymer 
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The dynamics of flexible polymer molecules are often assumed to be governed by hydrodynamics of 
the solvent. However there is considerable evidence that internal dissipation of a polymer contributes 
as well. Here we investigate the dynamics of a single chain in the absence of solvent to characterize 
the nature of this internal friction. We model the chains as freely hinged but with localized bond 
angles and 3-fold symmetric dihedral angles. We show that the damping is close but not identical 
to Kelvin damping, which depends on the first temporal and second spatial derivative of monomer 
position. With no internal potential between monomers, the magnitude of the damping is small for 
long wavelengths and weakly damped oscillatory time dependent behavior is seen for a large range of 
spatial modes. When the size of the internal potential is increased, such oscillations persist, but the 
damping becomes larger. However underdamped motion is present even with quite strong dihedral 
barriers for long enough wavelengths. 

PACS numbers: 



I. INTRODUCTION 

Polymer molecules in solution are usually modeled as- 
suming that damping is mediated by interaction with sur- 
rounding fluid, or other polymer chains. However it has 
been noted early on, that there is another source of dis- 
sipation, internally generated by the nonlinear dynamics 
of the polymer chain [iHl]- In a set of important ex- 
periments 3 extensional relaxation for different solvent 
viscosities was measured. By extrapolating the solvent 
viscosities to zero, an interesting residual dissipation was 
discovered that is internal in origin. A number of expla- 
nations for it have been proposed d, Q , but the origins 
of this effect are still not well understood. 

In this work, the internal friction of a polymer is ex- 
amined by considering it in isolation, that is, without 
any solvent. This not only sheds light on the case of 
polymer solutions, but also on more recent experimental 
techniques, discussed below, used to characterize poly- 
mers by first gasifying them without damaging their in- 
tegrity. 

Developments in mass spectrometry of long chain poly- 
mers, have become important in recent years in order 
to characterize large biological molecules such as pro- 
teins [Hi . The process involved in these experiments puts 
single chains into a vacuum. The study of polymer damp- 
ing in this environment would elucidate the understand- 
ing of internal dissipation but has yet to be studied exper- 
imentally. However recent work has studied the behavior 
of single chain molecules in a vacuum theoretically and by 
means of computer simulation The statistical prop- 

erties of single chains with no damping, and subject only 
Newton's laws was investigated. Both ideal chains (where 
only chain connectivity is included and no intra-chain 
interactions), and self-avoiding chains, were considered. 
In related work, the detailed chaotic properties of small 
chains with rigid links has also been observed for self- 
interacting chains with two body Lennard- Jones poten- 
tials . Through a comprehensive analysis of this prob- 



lem, they were able to show that the simulated chains 
were in excellent agreement with exact predictions from 
a microcanonical average, providing strong evidence that 
the dynamics are indeed ergodic. The chaotic motion 
derives from both the bond constraints and the convex 
nature of scattering. In another related work [l3| com- 
puter simulations for polymer molecules were performed 
to examine Lyapunov exponents. These were shown to 
be very small for systems close to a phase transition, such 
as the coil-globule transition. 

An ideal chain modeled with only linear springs con- 
necting adjacent monomers is integrable and shows no 
damping. Energy in every mode is conserved and cannot 
be exchanged with other modes. Dissipation results from 
the nonlinear coupling of modes to each other. This is 
closely related to the problem of nonlinear one dimen- 
sional chains that have been studied extensively using 
the Fermi, Pasta, and Ulam (FPU) model [lH. In this 
model, there is energy transfer between modes, however 
depending on initial conditions, the time it takes to lose 
correlation with its initial state can be very long . An 
ideal chain is a similar one dimensional system where the 
displacement of each particle is, in this case, much larger 
and is vectorial in nature. However we expect that the 
same equilibration problems persist. Therefore to ther- 
malize a polymer in a vacuum efficiently, the model cho- 
sen [7| uses rigid links so that the local motion would 
be highly chaotic. It also mimics the fact that the fluc- 
tuations in bond lengths are small compared to those 
involved with rotational degrees of freedom. On large 
length scales however, one might expect universal behav- 
ior, independent of the form of the potential connecting 
adjacent monomers. As will be seen below, this appears 
to be the case. 

The dynamics of such ideal chains are also closely re- 
lated to the study of the dynamics and energy flow in one 
dimensional chains, that govern its anomalous conductiv- 
ity [l3| • The dynamics in one dimension systems with 
momentum conservation show faster relaxation than in 
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higher dimensions due to Gahlean invariance. A con- 
sequence of momentum conservation in higher dimen- 
sions is responsible for long time tails in liquids that 
were predicted theoretically \lE\ and confirmed experi- 
mentally [13 . In the case of one dimensional nonlinear 
chains, a wide range of different models give the same 
heat conductivity exponent, for example the Sinai Pen- 
case model , the Random Collision Model [3] and 
FPU chains [18 1. The main difference between these one 
dimensional systems and ideal chains in a vacuum is that 
the former has a fractal dimension of 1, and the latter, of 
2. It is therefore hard to define local hydrodynamic vari- 
ables as a function of position for polymer chains. How- 
ever conservation of energy and momentum are common 
to both problems. 

In addition, in the case of an ideal chain in a vacuum, 
the equilibrium properties are effected by the conserva- 
tion of angular momentum and this can be analyzed 
exactly, showing that when the angular momentum is 
zero, the radius of gyration is significantly smaller than 
without that conservation law enforced. 

Initial results Q for an ideal chain in a vacuum showed 
that its dynamics are very different than those of a chain 
in solution. A freely hinged chain of N links, with con- 
stant link lengths was studied. The time auto-correlation 
function for position oscillates and is slowly damped. The 
damping time appears to scale as Trei oc ]\[(i-&5±.i5) ^ 
where L is the chain length. 

The reason for this behavior can be seen by under- 
standing the form of dissipation such chains should have. 
Because of Galilean invariance, the frictional force fd on 
a monomer cannot be proportional to its velocity, as this 
would imply that a uniformly translating chain would 
slow down. If we denote the position of a chain at ar- 
clength s as r(s), then the simplest term respecting this 
translational invariance is 



fd oc 



q3j. 

dtdh 



(1) 



This is the form of "Kelvin Damping" [19]. This theo- 
retical form was proposed by Maclnnes based on a cal- 
culation done by him for a model system It was 
later proposed as the origin of Cerf friction 0] where 
it was argued that this form was compatible with ex- 
periments characterizing internal friction Q. In Fourier 
space with wave-vector k conjugate to s, and oj conju- 
gate to t, this dissipation is proportional to iujk'^r{k, w). 
Longer wavelength modes are weakly damped, imply- 
ing underdamped motion for long enough wavelengths 
as we will shortly see. Using this damping along with 
linear forces between monomers, analogous to the Rouse 
model [i^l, one has 
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q2^ 



Kr"{s) is identical to that of the Rouse equation and de- 
scribes the net force experienced from neighboring seg- 
ments. For Gaussian chains, k = iksT, where T is the 
temperature For a linear chain, r'(s) = at the two 
ends and the chain can be expanded in terms of Fourier 
modes, similar to the analysis of the Rouse equation 



r(s) = ^ffc cos(fcs) 



with k = rnr/L, n = 1,2 

q2^ 



so that 
d 



(3) 



(4) 



When noise is omitted, this is of the form a damped 
harmonic oscillator with an effective mass of p, a spring 
constant oi K — nk"^ and friction of = Cfc^ , so that 



Mr ^- vr + Kr = 



(5) 



Here ^ is a noise term added to maintain ideal chain 
statistics, p is the mass per unit arclength. The term 



For underdamped long wavelength modes, the solution 
is proportional to exp(int) with flk = i-^k + i^k- The 
correlation function for a k mode can then be calculated 



{rk{0)rl{t)) = {\rkf)Re (exp(il]fci)^^fcM)- (6) 

Afc oc k^ (in the underdamped regime) so the damping 
time oc 1/fc^. The frequency of oscillation, w^, is propor- 
tional to k for small k . Therefore in this limit, there are 
many oscillations, 0{k~^) in a damping time. 

The model of a freely hinged chain is not realistic at 
a microscopic level. A molecule such as polyethylene has 
a strong orientational dependence as dihedral angles are 
varied. The chain will spend most of its time in min- 
ima and make transitions between these. The arguments 
above are general and do not give any indication of the 
prefactor for the dissipation term. This is expected to 
depend on the details of the potential. We will use nu- 
merical methods to see to what extent Eq. 2] is satisfied 
and to determine how the size of the prefactor C depends 
on the potential used. We will find, surprisingly, that un- 
derdamped motion is still present even with rather strong 
potentials. We first describe the numerical method used 
in this work and then give the simulation results in the 
following section. 



II. NUMERICAL METHOD 

This method is an extension of the simulation of rigid 
link systems developed by the author previously (2ll . [22| 
which considered the case of a highly overdamped sys- 
tems where inertia was negligible. Here we consider the 
general case of particles with mass m and damping 7. In 
the subsequent sections, we will take 7 = 0. 

The coordinates are denoted ri, r2, . . . , rjv, where iV is 
the number of masses in the system. The time derivative 
of these, that is the velocities, are denoted Vi, V2, . . . , v^y. 



3 



As in the earlier work, we will assume that each mass is 
being acted on by a force f^, which can be due to self- 
interaction or externally applied. It will be convenient 
to define the difference operator of adjacent coordinates 
or velocities A^r = Tj+i — and A^v = v^+i — v^. To 
keep the link lengths constant, we Introducing Lagrange 
multipliers ti,. . . ,tN~i which describe the tensions be- 
tween neighboring masses, we can write the equations of 
motion as: 

m-Vi+'fVi = ii(A^r) - ii_i(Aj_ir) + (7) 
r^ - V, (8) 

for 1 < i < iV. For a linear chain, one can define io = 
In =: to give the equations for the chain ends, that is 
for i = 1 and i — N. We will discuss ring chains below. 

We wish to evolve the r's and v's in time by finding the 
tension at every time step which allows us to iterate the 
above equations by some appropriate integration scheme. 
However this can be problematic due to the cumulation 
of errors, as we require that the magnitude of A^r remain 
very close to the step length, I. Define the error in this 
quantity as 

e, = |A,rp - P. (9) 

Without the problem of numerical error, the tensions 
would be determined by the condition = for i = 
1, . . . , N. By differentiating this equation twice with re- 
spect to time, a formula for the tensions can be obtained. 
However numerical error will cause bond lengths to stray 
from their initial values. Therefore we need to intro- 
duce feedback into the method so that non-zero will 
be pushed back towards zero. So instead we consider the 
equation 

Ae, + Bu + e; = (10) 

where A and B are constants causing a damping of errors 
with time, and whose values are determined to maximize 
computational efficiency. This can be rewritten as 

2A,r • A,v = -Ae, - 2BA,r • A,v - 2|A,vp (11) 

Applying the difference operator A^ to Eq. [7] and tak- 
ing the dot product with A^r gives 

mA^r • AiV -I- 7Air • A^v 

^U+iA^v ■ Ai+ir + ii-iA^r • A^^ir- 
2t,|A,r|2-f A,r- A,f (12) 

Putting this into Eq. [TT] gives 

<i+iA,r • A,+ir - 2U\1^,y\^ + U^i^v ■ Ai_ir 

771 

= -(-Ae,; - 2BA,r • A,v - 2|A,vn 

- A,r- A,f + 7A,r • A,v (13) 

The left hand side contain the tensions, which are un- 
knowns, but the rest of the variables and the right hand 



side are all known. These equations form a tridiagonal 
matrix equation which can be solved for the ti 's in a time 
0{N). 

After the tensions are determined, they are used in the 
right hand side of Eq. [T] These 2A^ first order differen- 
tial equation can be integrated by a variety of methods. 
This paper uses fourth order Runge Kutta to do this. 
Excluding the operations involving the computation of 
the forces f , this algorithm runs in 0{N) operations per 
time step. For short range potentials, such as are used 
here, the forces also require 0{N) operations, meaning 
the operations per time step are 0{N). This scales the 
same way with N as non-constrained simulations such as 
molecular dynamics with variable bond lengths. 

In general, this algorithm can be easily extended to sys- 
tems with any set of connections between masses, such 
as branched topology. For the simplest variant, the ring 
chain, we modify the problem by introducing fictitious 
particles r^v+i = ri, and Fq = rjy. This places an ad- 
ditional link between monomer 1 and N, giving rise to 
an additional tension In- To keep the form of the equa- 
tions the same, it is convenient to introduce tg = In- 
In this case Eq. [13] becomes a cyclic tridiagonal system 
of equations that can be easily solved as well. However 
for more general connections, the order of the number of 
operations will be greater than N. 

Typically, the step size used in this work for a Runge 
Kutta iteration was .01. Therefore this algorithm pro- 
vides an efficient method for investigating dynamics of a 
polymer in a vacuum. 



III. SIMULATIONS 

We first consider the case where there are no potentials 
but only the freely hinged constraint. In this work, we 
consider a step length / = 1 and at temperature ksT = 1. 
We fit the correlation functions for numerical data for 
iV = 64 using Eq. (5] The first two mode are shown in 
Fig. [T] As can be seen, the fit to the data is excellent. 
This was averaged over 37578 runs, so that the error 
bars on the data are negligible. Clearly for this range of 
parameters, the decay of a single fc-mode is well described 
by a damped harmonic oscillator. All data for correlation 
functions shown are normalized to unity at < = 0. This 
is because the amplitude of the modes is proportional to 
1/fc^, making it hard to discern the data without this 
rescaling. 

It is important to note that in general, this kind of cor- 
relation function can not be described by a single mode, 
and that the answer is expected to be the sum over a 
large number of modes. In fact, for strong enough dihe- 
dral potentials more modes need to be included, however 
we shall see that for quite strong potentials, a single mode 
fits the simulation data quite well. 

The parameters for the oscillator can be fit as a func- 
tion of k. flk is first determined and then the correspond- 
ing parameters in Eq. [Sj M/K and v/K are calculated. 
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FIG. 1: (Color Online) Correlation functions for a freely 
hinged chain in a vacuum for A'^ = 64 for the lowest two 
modes, shown by the solid triangles. The line going through 
them are fits using Eq. [B] 



In Fig. [21 M and v were fitted for different mode num- 
bers fc — Tin/N with N = 256 (the step length has been 
set to unity). The fluctuations give a measure of the 
uncertainty in this data as was checked by fitting with 
independent data. 
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FIG. 2: (Color Online) The parameters for fits of the time 
auto-correlation function, Eq. [6] to a damped harmonic oscil- 
lator are shown for a freely hinged chain with A'^ = 256. The 
upper line represents the mass k^M/K of each mode, and the 
lower curve represents that damping v/ K. 



Note that according to Eq. [51 both curves should be 
constant as a function of the mode number. However a 
clear variation is seen in each. The effective mass of the 
mode increases with increasing fc, and the effective fric- 
tion coefhcient is not oc fc^, but is slightly less. The data 



could be fit equally well, with a slightly smaller exponent 
so that v (X kP with p < 2, or as a logarithmic correc- 
tion, k'^ / log(fc). In either case, the assumption of a local 
term for the friction is apparently not completely cor- 
rect. There is a subtle transfer of energy between modes 
leading to non-locality. The same is apparently true for 
the mass, and no doubt these two effects are related. The 
goal of the present work is not to explain these effects the- 
oretically, but to investigate how the friction coefficient 
is altered by local potentials along the chain. Therefore 
we will now turn our attention to this problem. 

The oscillations seen are a result of the weak form of 
damping at low wave-number. Do these persist if the 
chain is no longer freely hinged? To answer this, we 
first examine a model where the distribution of bond, or 
valency angles, is weighted around a set of values centered 
at a particular angle Oy ■ To restrict configurations in this 
manner, a potential is constructed between the nearest 
neighbors of monomer i as follows: 

f/.ai = ^(|r.+i-r,_ip-r2)2. (14) 

Because the distance between nearest neighbors is fixed, 
for large Ub, this will limit configurations as just de- 
scribed. To is chosen to give the value of 9v desired, 
in this case 9v — 104°. This value is somewhat arbi- 
trary and varies according to the chemical structure of 
the polymer. As a result, the dynamics will be restricted 
as well. Fig. [3] shows the time auto-correlation function 
for Ub = 10, for the first three modes. In this case we 
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FIG. 3: (Color Online) Time auto-correlation functions 
where bond angle is now restricted for A'' — 64. Data for 
the first three modes is shown by the solid triangles. The 
lines going through it is a fit using Eq. |6| 



see again that the fit to Eq. [6| is excellent. The long 
wavelength modes are still quite underdamped. A fit of 
the harmonic oscillator parameters is shown in Fig. [31 
It is similar to the freely hinged case. Fig. [1] suggesting 
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that the variation of the drag coefScient with k may be 
understandable by some general mechanism. lo 
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FIG. 4: (Color Online) The parameters for fits of the time 
auto-correlation function, Eq. [S] to a damped harmonic os- 
cillator are shown for the same data as in Fig. [S] The upper 
curve is a measure of the mass k^M/K, and the lower curve 
is measure of the damping v/K. 

Usually as a link is rotated, there are three energy 
minima as a function of the dihedral angle. To make 
this model more realistic, it is necessary to add such a 
dihedral potential. As a function of the dihedral angle 
9, between two adjacent monomers, the potential energy 
U ~ Vd cos(36'). This can be expressed in terms of cos^ 
using the usual trigonometric identity. If the bond angle 
is fixed, cos6' can be expressed in terms of dot prod- 
ucts which is computationally advantageous. Because 
the bond angle is almost constant, we still use the same 
dot product formula in the simulation. This also breaks 
symmetry between the three minima moving the trans 
state slightly relative to the two gauche states. This is 
also seen in real data such as for polyethylene [1^. Fig. 
[5] shows the distribution of bond angles, p{0), for two 
separate dihedral potentials, with = 3 and = 4, 
for Ub = 10 and N = 64. As the amplitude of the 
dihedral potential, Vd increases, the ratio of the maxi- 
mum for the distribution, pmax to the minimum pmin 
increases. In these simulations the energy scale was cho- 
sen so that the temperature is 1 and Fig. [6] shows that 
^og{pmax / Pmin) ~ 2Vd, It is reduccd from approximately 
4Vd because of coupling to other degrees of freedom such 
as bond angle. For a real temperature of SSO-fiT, compari- 
son with earlier modeling [23] of polyethylene, yields that 
Vd ~ 2.5 (as we are using units where T = 1.) At lower 
real temperatures, the potential will be correspondingly 
higher. 

Fig. [7] shows the correlation function fit to the damped 
harmonic oscillator for Vd = 2 (a), and Vd — 3 (b). Here 
iV = 64 and Ub = 10. In both cases, clear oscillations 
can be seen for the lowest modes. 



FIG. 5: (Color Online) The distribution of dihedral bond 
angles for two different heights of the dihedral potential, Vd = 
3 (circles) and Vd = 4: (triangles). Here Ub = 10 and — 64. 
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FIG. 6: (Color Online) The natural logarithm of the dihe- 
dral angle distribution maximum to minimum as a function 
of potential amplitude Vd with Ub ~ 10. 



It is of interest to know how the frequency dependent 
damping depends on the height of the dihedral poten- 
tial. Fig. |S] shows the results of fitting different modes 
to the damped harmonic oscillator, taking into account 
the possibility of both underdamped and overdamped 
motion. As with Fig. [2l the damping is divided by 
K cc k^. The non-constant nature of the results show 
that for quite substantial Vd, there is a similar deviation 
from the expected Kelvin damping as was seen in Figs. [2] 
and[31 Note that even if higher modes are overdamped we 
expect the for lower enough fc, and long enough chains, 
the motion will become underdamped. This is because 
of the /c-dependent form of the damping v, which goes to 
zero for as fc ^> 0. 

As can be seen from Fig. HI for > 1, barriers to 
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FIG. 7: (Color Online) Time auto-correlation functions 
where bond angle and a dihedral potential for TV — 64, 
Ub = 10. Data for the first three modes is shown by the 
solid triangles. The lines going through it is a fit using Eq. [S] 
(a) with Vd = 2 and (b) with Vd = 3. 



rotation are substantial. It would then first appear that 
motion would be much more like that of a lattice model, 
where a monomer stays in a potential minimum for a long 
time and flips to another state. This would suggest that 
the motion would be heavily overdamped, with correla- 
tion functions like those of the Rouse model. However 
we have seen that even in these cases, the motion for 
long wavelengths, is underdamped. For strong enough 
potentials such as Vd = 4, the motion even at long wave- 
lengths and N = 64, ceases to be underdamped. The 
relaxation becomes very slow and does not fit well to a 
simple harmonic oscillator. Many more frequency modes 
must be considered in this case. However as we argued 



above, for long enough chain length, we expect to see 
underdamped motion for low k. When the value of Ub 
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FIG. 8: (Color Online) The damping as a function of mode 
number for a range of dihedral potentials with A'^ = 64. The 
bottom curve shown by the triangles is Vd = 0, the next 
highest, Vd = 1 (squares), then Vd ~ 2 (diamonds), and Vd = 
3 (crosses). 



is increased, this increases \og{pmax / Pmin) and will make 
the barrier height larger. Still for N = 64, when the dihe- 
dral constant Ub is raised to J7b = 40, oscillations persist 
when = 2. 

IV. CONCLUSIONS 

These results are of interest for two reasons. First, it 
sheds like on the nature of internal friction of polymer 
chains in solution. This work characterizes the nature 
of this dissipation. It is quite similar to that of Kelvin 
friction, Eq. [TJ However it is not identical as the damp- 
ing is larger for small wavenumber. This departure is at 
present unexplained and is probably quite nontrivial to 
understand as it is related to the dynamics of nonlinear 
one dimensional chains [ill, US] • 

Second, it is important in the understanding the inter- 
nal dynamics of polymers in a vacuum of which there is 
at present little experimental or theoretical understand- 
ing. In reality polymers in this situation will be charged 
and have van der Waals interactions. An investigation 
of these have shown that their effects are very impor- 
tant 0. However as with the understanding of ideal 
chains in polymer solutions, it is important to under- 
stand an ideal chain in a vacuum as the starting point 
for further analysis. 

The author thanks Onuttom Narayan for useful dis- 
cussions. 
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